Network Fitting
Emily Pollock
30 December, 2020
1 Session Info
sessioninfo::package_info()## package * version date lib source
## ape 5.3 2019-03-17 [1] CRAN (R 3.6.0)
## assertthat 0.2.1 2019-03-21 [1] CRAN (R 3.6.0)
## backports 1.1.9 2020-08-24 [1] CRAN (R 3.6.2)
## cli 2.0.2 2020-02-28 [1] CRAN (R 3.6.0)
## coda 0.19-3 2019-07-05 [1] CRAN (R 3.6.0)
## codetools 0.2-16 2018-12-24 [1] CRAN (R 3.6.1)
## colorspace 1.4-1 2019-03-18 [1] CRAN (R 3.6.0)
## crayon 1.3.4 2017-09-16 [1] CRAN (R 3.6.0)
## DEoptimR 1.0-8 2016-11-19 [1] CRAN (R 3.6.0)
## deSolve * 1.27.1 2020-01-02 [1] CRAN (R 3.6.0)
## digest 0.6.25 2020-02-23 [1] CRAN (R 3.6.0)
## doParallel 1.0.15 2019-08-02 [1] CRAN (R 3.6.0)
## dplyr 1.0.2 2020-08-18 [1] CRAN (R 3.6.2)
## egonet * 0.0.0.9000 2020-02-11 [1] Github (statnet/egonet@70d4aa2)
## ellipsis 0.3.1 2020-05-15 [1] CRAN (R 3.6.2)
## EpiModel * 2.0.3 2020-11-09 [1] CRAN (R 3.6.2)
## ergm * 3.11.0 2020-10-14 [1] CRAN (R 3.6.2)
## ergm.ego * 0.6.1 2020-11-19 [1] CRAN (R 3.6.2)
## evaluate 0.14 2019-05-28 [1] CRAN (R 3.6.0)
## fansi 0.4.1 2020-01-08 [1] CRAN (R 3.6.0)
## foreach 1.4.7 2019-07-27 [1] CRAN (R 3.6.0)
## generics 0.0.2 2018-11-29 [1] CRAN (R 3.6.0)
## ggplot2 3.3.2 2020-06-19 [1] CRAN (R 3.6.2)
## glue 1.4.1 2020-05-13 [1] CRAN (R 3.6.2)
## gtable 0.3.0 2019-03-25 [1] CRAN (R 3.6.0)
## here * 0.1 2017-05-28 [1] CRAN (R 3.6.0)
## htmltools 0.5.0 2020-06-16 [1] CRAN (R 3.6.2)
## iterators 1.0.12 2019-07-26 [1] CRAN (R 3.6.0)
## knitr 1.29 2020-06-23 [1] CRAN (R 3.6.2)
## lattice 0.20-38 2018-11-04 [1] CRAN (R 3.6.1)
## lazyeval 0.2.2 2019-03-15 [1] CRAN (R 3.6.0)
## lifecycle 0.2.0 2020-03-06 [1] CRAN (R 3.6.0)
## lpSolve 5.6.13.3 2019-08-19 [1] CRAN (R 3.6.0)
## magrittr 1.5 2014-11-22 [1] CRAN (R 3.6.0)
## MASS 7.3-51.4 2019-03-31 [1] CRAN (R 3.6.1)
## Matrix 1.2-17 2019-03-22 [1] CRAN (R 3.6.1)
## munsell 0.5.0 2018-06-12 [1] CRAN (R 3.6.0)
## network * 1.16.1 2020-10-07 [1] CRAN (R 3.6.2)
## networkDynamic * 0.10.0 2019-04-05 [1] CRAN (R 3.6.0)
## nlme 3.1-140 2019-05-12 [1] CRAN (R 3.6.1)
## pillar 1.4.6 2020-07-10 [1] CRAN (R 3.6.2)
## pkgconfig 2.0.3 2019-09-22 [1] CRAN (R 3.6.0)
## purrr 0.3.4 2020-04-17 [1] CRAN (R 3.6.2)
## R6 2.4.1 2019-11-12 [1] CRAN (R 3.6.0)
## RColorBrewer 1.1-2 2014-12-07 [1] CRAN (R 3.6.0)
## Rcpp 1.0.5 2020-07-06 [1] CRAN (R 3.6.2)
## rlang 0.4.7 2020-07-09 [1] CRAN (R 3.6.2)
## rle 0.9.2 2020-09-25 [1] CRAN (R 3.6.2)
## rmarkdown 2.3 2020-06-18 [1] CRAN (R 3.6.2)
## robustbase 0.93-5 2019-05-12 [1] CRAN (R 3.6.0)
## rprojroot 1.3-2 2018-01-03 [1] CRAN (R 3.6.0)
## scales 1.1.1 2020-05-11 [1] CRAN (R 3.6.2)
## sessioninfo 1.1.1 2018-11-05 [1] CRAN (R 3.6.0)
## statnet.common 4.4.1 2020-10-03 [1] CRAN (R 3.6.2)
## stringi 1.4.6 2020-02-17 [1] CRAN (R 3.6.0)
## stringr 1.4.0 2019-02-10 [1] CRAN (R 3.6.0)
## tergm * 3.6.1 2019-06-12 [1] CRAN (R 3.6.0)
## tergmLite * 2.2.1 2020-07-22 [1] CRAN (R 3.6.2)
## tibble 3.0.3 2020-07-10 [1] CRAN (R 3.6.2)
## tidyselect 1.1.0 2020-05-11 [1] CRAN (R 3.6.2)
## trust 0.1-8 2020-01-10 [1] CRAN (R 3.6.0)
## vctrs 0.3.2 2020-07-15 [1] CRAN (R 3.6.2)
## withr 2.2.0 2020-04-20 [1] CRAN (R 3.6.2)
## xfun 0.16 2020-07-24 [1] CRAN (R 3.6.2)
## yaml 2.2.1 2020-02-01 [1] CRAN (R 3.6.0)
##
## [1] /Library/Frameworks/R.framework/Versions/3.6/Resources/library2 General Controls
need to update duration targets w/ estimated means from parametric model
but then will also need to use CKs correction in-simiulation find undelrying relationship length in presense of mortality
np <- detectCores()
# Set parallel.type
if (np==0) ptype=NULL else ptype='PSOCK'
dat <- setup_params()
ppopsize <- dat$ppopsize
MCMCinterval <- dat$mcmcInterval
MCMCburnin <- dat$mcmcBurnin
time.step <- dat$time.step
age.width <- dat$age.width
dissolution <- dat$dissolution
#mRate <- dat$mRate
duration.mar <- dat$duration.marriage
duration.cohab <- dat$duration.cohab
duration.other <- dat$duration.other
DissolutionMar <- dissolution_coefs(dissolution=dissolution,
duration=duration.mar)
DissolutionCohab <- dissolution_coefs(dissolution=dissolution,
duration=duration.cohab)
DissolutionOther <- dissolution_coefs(dissolution=dissolution,
duration=duration.other)3 Some Adjustments
- make age/sex attribute (i.e. 15F, 23M)
- sample down to 10,000
ssize <- 10000
# sample down to 10,000
egodat_marriage <- sample(egodat_marriage, ssize)
# pull sampled egos, attribute updates, and update other alter lists (for consistent nodeset)
egos <- egodat_marriage$egos
egodat_cohab$alters <- egodat_cohab$alters[egodat_cohab$alters$ego %in% egos$ego,]
egodat_other$alters <- egodat_other$alters[egodat_other$alters$ego %in% egos$ego,]
egodat_once$alters <- egodat_once$alters[egodat_once$alters$ego %in% egos$ego,]
egos$ageF <- ifelse(egos$male==0, egos$age, 0)
egos$ageM <- ifelse(egos$male==1, egos$age, 0)
egodat_marriage$alters$ageF <- ifelse(egodat_marriage$alters$male==0, egodat_marriage$alters$age, 0)
egodat_marriage$alters$ageM <- ifelse(egodat_marriage$alters$male==1, egodat_marriage$alters$age, 0)
egodat_cohab$alters$ageF <- ifelse(egodat_cohab$alters$male==0, egodat_cohab$alters$age, 0)
egodat_cohab$alters$ageM <- ifelse(egodat_cohab$alters$male==1, egodat_cohab$alters$age, 0)
egodat_other$alters$ageF <- ifelse(egodat_other$alters$male==0, egodat_other$alters$age, 0)
egodat_other$alters$ageM <- ifelse(egodat_other$alters$male==1, egodat_other$alters$age, 0)
egodat_once$alters$ageF <- ifelse(egodat_once$alters$male==0, egodat_once$alters$age, 0)
egodat_once$alters$ageM <- ifelse(egodat_once$alters$male==1, egodat_once$alters$age, 0)
egodat_marriage$egos <- egos
egodat_cohab$egos <- egos
egodat_other$egos <- egos
egodat_once$egos <- egos
# grab binary deg dist for plots below, x2 for scaline to 20000 nodes
m <- with(egodat_marriage$egos[egodat_marriage$egos$male==1,], table(age, deg.mar))
m <- m[,2]*(ppopsize/ssize)
f <- with(egodat_marriage$egos[egodat_marriage$egos$male==0,], table(age, deg.mar))
f <- f[,2]*(ppopsize/ssize)
mc <- with(egodat_cohab$egos[egodat_cohab$egos$male==1,], table(age, deg.cohab))
mc <- mc[,2]*(ppopsize/ssize)
fc <- with(egodat_cohab$egos[egodat_cohab$egos$male==0,], table(age, deg.cohab))
fc <- fc[,2]*(ppopsize/ssize)
mo <- with(egodat_other$egos[egodat_other$egos$male==1,], table(age, deg.other))
mo <- mo[,2]*(ppopsize/ssize)
fo <- with(egodat_other$egos[egodat_other$egos$male==0,], table(age, deg.other))
fo <- fo[,2]*(ppopsize/ssize)4 Marriage Network
4.1 Ergm.Ego Estimation & Setup
fit.marriage <- ergm.ego(egodat_marriage ~ edges +
nodecov(~ageF) +
nodecov(~ageF^2) +
nodecov(~ageF<30) +
nodecov(~ageF>30) +
nodecov(~ageM) +
nodecov(~ageM^2) +
nodecov(~ageM<30) +
nodecov(~ageM>30) +
absdiff(~sqrtage + dat$shift.mar*(male==0)) +
nodefactor("deg.other.binary") +
offset(nodematch("male", diff = FALSE)) +
offset(nodefactor("olderpartnerM")) +
#offset(nodefactor("debuted", levels=1)) +
offset("concurrent"),
offset.coef = c(-Inf, -Inf, -Inf), #offset.coef = c(-Inf, -Inf, -Inf, -Inf),
control =
control.ergm.ego(ppop.wt = "sample", ppopsize = ppopsize,
ergm.control =
control.ergm(MCMLE.maxit = 100,
MCMC.interval = MCMCinterval,
MCMC.burnin = MCMCburnin,
parallel = np,
parallel.type = ptype)))
#nodefactor(~ageF, levels=5:9) +
#nodefactor(~ageM, levels=7:10) +
saveRDS(fit.marriage, here("fits", "fit.marriage.rds"))fit.marriage <- readRDS(here("fits", "fit.marriage.rds"))
mpop <- environment(fit.marriage$ergm.formula)$popnw
marriage.netest <- ego.netest(fit.marriage, DissolutionMar)
summary(fit.marriage)## Call:
## ergm(formula = ergm.formula, offset.coef = ergm.offset.coef,
## target.stats = m, eval.loglik = FALSE, control = control$ergm.control)
##
## Iterations: 3 out of 100
##
## Monte Carlo MLE Results:
## Estimate Std. Error MCMC % z value
## offset(netsize.adj) -9.903488 0.000000 0 -Inf
## edges -18.415798 2.313703 0 -7.959
## nodecov.ageF 0.658567 0.103593 0 6.357
## nodecov.ageF^2 -0.008985 0.001596 0 -5.630
## nodecov.ageF<30 -0.302788 0.503616 0 -0.601
## nodecov.ageF>30 -0.696689 0.513842 0 -1.356
## nodecov.ageM 0.581784 0.113492 0 5.126
## nodecov.ageM^2 -0.007603 0.001686 0 -4.509
## nodecov.ageM<30 0.381250 0.295776 0 1.289
## nodecov.ageM>30 0.260988 0.294664 0 0.886
## absdiff.sqrtage+dat$shift.mar*(male==0) -3.700792 0.141901 0 -26.080
## nodefactor.deg.other.binary.1 -4.394601 0.326862 0 -13.445
## offset(nodematch.male) -Inf 0.000000 0 -Inf
## offset(nodefactor.olderpartnerM.1) -Inf 0.000000 0 -Inf
## offset(concurrent) -Inf 0.000000 0 -Inf
## Pr(>|z|)
## offset(netsize.adj) <1e-04 ***
## edges <1e-04 ***
## nodecov.ageF <1e-04 ***
## nodecov.ageF^2 <1e-04 ***
## nodecov.ageF<30 0.548
## nodecov.ageF>30 0.175
## nodecov.ageM <1e-04 ***
## nodecov.ageM^2 <1e-04 ***
## nodecov.ageM<30 0.197
## nodecov.ageM>30 0.376
## absdiff.sqrtage+dat$shift.mar*(male==0) <1e-04 ***
## nodefactor.deg.other.binary.1 <1e-04 ***
## offset(nodematch.male) <1e-04 ***
## offset(nodefactor.olderpartnerM.1) <1e-04 ***
## offset(concurrent) <1e-04 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## The following terms are fixed by offset and are not estimated:
## offset(netsize.adj) offset(nodematch.male) offset(nodefactor.olderpartnerM.1) offset(concurrent)4.2 MCMC Diagnostics / GOF
mcmc.diagnostics(fit.marriage)
gof <- gof(fit.marriage, GOF = "model")
saveRDS(gof, here("fits", "gof.m.rds"))gof <- readRDS(here("fits", "gof.m.rds"))
plot(gof)
4.3 Static Diagnostics
static.mar <- netdx(marriage.netest, nsims = 500, dynamic = FALSE)
saveRDS(static.mar, here("fits", "static.m.rds"))static.mar <- readRDS(here("fits", "static.m.rds"))
static.mar## EpiModel Network Diagnostics
## =======================
## Diagnostic Method: Static
## Simulations: 500
##
## Formation Diagnostics
## -----------------------
## Target Sim Mean Pct Diff
## edges 3478.196 3474.476 -0.107
## nodecov.ageF 115489.827 115372.034 -0.102
## nodecov.ageF^2 3960149.254 3956455.334 -0.093
## nodecov.ageF<30 4512.671 4506.556 -0.136
## nodecov.ageF>30 2282.498 2280.464 -0.089
## nodecov.ageM 120933.300 120815.776 -0.097
## nodecov.ageM^2 4326955.505 4322977.556 -0.092
## nodecov.ageM<30 4216.771 4211.118 -0.134
## nodecov.ageM>30 2562.346 2561.888 -0.018
## absdiff.sqrtage+dat$shift.mar*(male==0) 819.197 820.059 0.105
## nodefactor.deg.other.binary.1 24.909 24.984 0.303
## offset(nodematch.male) NA 0.000 NA
## offset(nodefactor.olderpartnerM.1) NA 0.000 NA
## offset(concurrent) NA 0.000 NA
## Sim SD
## edges 29.796
## nodecov.ageF 969.749
## nodecov.ageF^2 34934.018
## nodecov.ageF<30 46.456
## nodecov.ageF>30 23.726
## nodecov.ageM 1001.210
## nodecov.ageM^2 36557.584
## nodecov.ageM<30 45.093
## nodecov.ageM>30 25.112
## absdiff.sqrtage+dat$shift.mar*(male==0) 15.012
## nodefactor.deg.other.binary.1 4.965
## offset(nodematch.male) 0.000
## offset(nodefactor.olderpartnerM.1) 0.000
## offset(concurrent) 0.000plot(static.mar)
4.4 Dynamic Diagnostics
dynamic.marriage <- netdx(marriage.netest, nsims = 5, nsteps = 1000, ncores = np,
set.control.ergm = control.simulate.ergm(MCMC.burnin=1e5,
MCMC.interval=1e5),
set.control.stergm = control.simulate.network(MCMC.burnin.min=1e5,
MCMC.burnin.max=1e5))
saveRDS(dynamic.marriage, here("fits", "dynamic.m.rds"))dynamic.marriage <- readRDS(here("fits", "dynamic.m.rds"))
plot(dynamic.marriage)
plot(dynamic.marriage)
plot(dynamic.marriage, type = "duration")
plot(dynamic.marriage, type = "dissolution")
4.5 Checking Mean Deg Dist by Age/Sex and Age/Sex Mixing are close to data
n <- dynamic.marriage$nw
plot(summary(n ~ nodefactor(~ageM)), type="l", main="Males")
lines(m, pch=19, col="blue")
legend("topright", legend = c("Empirical", "Model"), col=c("blue", "black"), lwd=1)
plot(summary(n ~ nodefactor(~ageF)), type="l", main="Females")
lines(f, col="blue")
legend("topright", legend = c("Empirical", "Model"), col=c("blue", "black"), lwd=1)
# age/sex mixing
aF <- n %v% "ageF"
aM <- n %v% "ageM"
el <- as.edgelist(n)
fem1 <- aF[el[,1]]
fem2 <- aF[el[,2]]
male1 <- aM[el[,1]]
male2 <- aM[el[,2]]
allFs <- NULL
allMs <- NULL
for (i in 1:length(fem1)){
if (fem1[i] != 0) {allFs[i] <- fem1[i]}
else {allFs[i] <- fem2[i]}
if (male1[i] != 0) {allMs[i] <- male1[i]}
else {allMs[i] <- male2[i]}
}
agesel <- as.data.frame(cbind(allFs, allMs))
sqrtages <- sqrt(agesel)
sqrtages$diff <- sqrtages[,1] - sqrtages[,2]
mean(sqrtages$diff)## [1] -0.13198955 Cohab Network
5.1 Ergm.Ego Estimation & Setup
x <- control.ergm.ego(ppopsize = 1, ppop.wt = "sample",
ergm.control = control.ergm(MCMLE.maxit = 200,
MCMC.interval = MCMCinterval,
MCMC.burnin = MCMCburnin))
x$ppopsize <- mpop
fit.cohab <- ergm.ego(egodat_cohab ~ edges +
nodecov(~ageF) +
nodecov(~ageF^2) +
nodecov(~ageM) +
nodecov(~ageM^2) +
nodefactor(~ageF, levels=5:9) +
nodefactor(~ageM, levels=7:10) +
absdiff(~sqrtage + dat$shift.cohab*(male==0)) +
nodefactor("deg.other.binary") +
offset(nodematch("male", diff = FALSE)) +
offset(nodefactor("olderpartnerC")) +
#offset(nodefactor("debuted", levels=1)) +
offset("concurrent"),
offset.coef = c(-Inf, -Inf, -Inf), #offset.coef = c(-Inf, -Inf, -Inf, -Inf),
control = x)
saveRDS(fit.cohab, here("fits", "fit.cohab.rds"))fit.cohab <- readRDS(here("fits", "fit.cohab.rds"))
cohab.netest <- ego.netest(fit.cohab, DissolutionCohab)
summary(fit.cohab)## Call:
## ergm(formula = ergm.formula, offset.coef = ergm.offset.coef,
## target.stats = m, eval.loglik = FALSE, control = control$ergm.control)
##
## Iterations: 2 out of 200
##
## Monte Carlo MLE Results:
## Estimate Std. Error MCMC % z value
## offset(netsize.adj) -9.9034876 0.0000000 0 -Inf
## edges -8.7874426 0.9013938 0 -9.749
## nodecov.ageF 0.4025384 0.0667472 0 6.031
## nodecov.ageF^2 -0.0077624 0.0012306 0 -6.308
## nodecov.ageM 0.3009793 0.0578495 0 5.203
## nodecov.ageM^2 -0.0046020 0.0009947 0 -4.627
## nodefactor.ageF.18 0.1432191 0.3019090 0 0.474
## nodefactor.ageF.19 0.6184190 0.2804912 0 2.205
## nodefactor.ageF.20 1.5189275 0.3457406 0 4.393
## nodefactor.ageF.21 0.8620661 0.2534431 0 3.401
## nodefactor.ageF.22 0.6722736 0.2334703 0 2.879
## nodefactor.ageM.20 0.5421379 0.3020196 0 1.795
## nodefactor.ageM.21 0.3026223 0.2997873 0 1.009
## nodefactor.ageM.22 0.7057102 0.2506985 0 2.815
## nodefactor.ageM.23 0.8937318 0.2602071 0 3.435
## absdiff.sqrtage+dat$shift.cohab*(male==0) -2.7147838 0.1335418 0 -20.329
## nodefactor.deg.other.binary.1 -2.4469524 0.3307137 0 -7.399
## offset(nodematch.male) -Inf 0.0000000 0 -Inf
## offset(nodefactor.olderpartnerC.1) -Inf 0.0000000 0 -Inf
## offset(concurrent) -Inf 0.0000000 0 -Inf
## Pr(>|z|)
## offset(netsize.adj) < 1e-04 ***
## edges < 1e-04 ***
## nodecov.ageF < 1e-04 ***
## nodecov.ageF^2 < 1e-04 ***
## nodecov.ageM < 1e-04 ***
## nodecov.ageM^2 < 1e-04 ***
## nodefactor.ageF.18 0.635230
## nodefactor.ageF.19 0.027470 *
## nodefactor.ageF.20 < 1e-04 ***
## nodefactor.ageF.21 0.000670 ***
## nodefactor.ageF.22 0.003983 **
## nodefactor.ageM.20 0.072647 .
## nodefactor.ageM.21 0.312756
## nodefactor.ageM.22 0.004878 **
## nodefactor.ageM.23 0.000593 ***
## absdiff.sqrtage+dat$shift.cohab*(male==0) < 1e-04 ***
## nodefactor.deg.other.binary.1 < 1e-04 ***
## offset(nodematch.male) < 1e-04 ***
## offset(nodefactor.olderpartnerC.1) < 1e-04 ***
## offset(concurrent) < 1e-04 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## The following terms are fixed by offset and are not estimated:
## offset(netsize.adj) offset(nodematch.male) offset(nodefactor.olderpartnerC.1) offset(concurrent)5.2 MCMC Diagnostics / GOF
mcmc.diagnostics(fit.cohab)
gof <- gof(fit.cohab, GOF = "model")
saveRDS(gof, here("fits", "gof.c.rds"))gof <- readRDS(here("fits", "gof.c.rds"))
plot(gof)
5.3 Static Diagnostics
static.cohab <- netdx(cohab.netest, nsims = 500, dynamic = FALSE)
saveRDS(static.cohab, here("fits", "static.c.rds"))static.cohab <- readRDS(here("fits", "static.c.rds"))
static.cohab## EpiModel Network Diagnostics
## =======================
## Diagnostic Method: Static
## Simulations: 500
##
## Formation Diagnostics
## -----------------------
## Target Sim Mean Pct Diff
## edges 1334.405 1333.460 -0.071
## nodecov.ageF 36428.502 36382.962 -0.125
## nodecov.ageF^2 1048047.957 1046062.370 -0.189
## nodecov.ageM 39792.620 39756.850 -0.090
## nodecov.ageM^2 1247690.001 1246314.762 -0.110
## nodefactor.ageF.18 26.033 25.994 -0.149
## nodefactor.ageF.19 58.637 58.818 0.309
## nodefactor.ageF.20 81.757 81.176 -0.711
## nodefactor.ageF.21 88.739 89.436 0.786
## nodefactor.ageF.22 83.736 84.506 0.919
## nodefactor.ageM.20 40.588 40.286 -0.745
## nodefactor.ageM.21 36.329 36.610 0.774
## nodefactor.ageM.22 89.100 89.390 0.325
## nodefactor.ageM.23 76.466 76.920 0.594
## absdiff.sqrtage+dat$shift.cohab*(male==0) 418.638 417.453 -0.283
## nodefactor.deg.other.binary.1 60.629 60.448 -0.298
## offset(nodematch.male) NA 0.000 NA
## offset(nodefactor.olderpartnerC.1) NA 0.000 NA
## offset(concurrent) NA 0.000 NA
## Sim SD
## edges 30.870
## nodecov.ageF 872.239
## nodecov.ageF^2 27193.786
## nodecov.ageM 948.242
## nodecov.ageM^2 31918.371
## nodefactor.ageF.18 4.686
## nodefactor.ageF.19 6.923
## nodefactor.ageF.20 6.859
## nodefactor.ageF.21 8.038
## nodefactor.ageF.22 7.621
## nodefactor.ageM.20 5.809
## nodefactor.ageM.21 5.477
## nodefactor.ageM.22 8.211
## nodefactor.ageM.23 7.485
## absdiff.sqrtage+dat$shift.cohab*(male==0) 14.249
## nodefactor.deg.other.binary.1 7.512
## offset(nodematch.male) 0.000
## offset(nodefactor.olderpartnerC.1) 0.000
## offset(concurrent) 0.000plot(static.cohab)
5.4 Dynamic Diagnostics
dynamic.cohab <- netdx(cohab.netest, nsims = 5, nsteps = 1000, ncores = np,
set.control.ergm = control.simulate.ergm(MCMC.burnin=1e5,
MCMC.interval=1e5),
set.control.stergm = control.simulate.network(MCMC.burnin.min=1e5,
MCMC.burnin.max=1e5))
saveRDS(dynamic.cohab, here("fits", "dynamic.c.rds"))dynamic.cohab <- readRDS(here("fits", "dynamic.c.rds"))
plot(dynamic.cohab)
plot(dynamic.cohab)
plot(dynamic.cohab, type = "duration")
plot(dynamic.cohab, type = "dissolution")
5.5 Checking Mean Deg Dist by Age/Sex and Age/Sex Mixing are close to data
n <- dynamic.cohab$nw
plot(summary(n ~ nodefactor(~ageM)), type="l", main="Males")
lines(mc, pch=19, col="blue")
legend("topright", legend = c("Empirical", "Model"), col=c("blue", "black"), lwd=1)
plot(summary(n ~ nodefactor(~ageF)), type="l", main="Females")
lines(fc, col="blue")
legend("topright", legend = c("Empirical", "Model"), col=c("blue", "black"), lwd=1)
# age/sex mixing
aF <- n %v% "ageF"
aM <- n %v% "ageM"
el <- as.edgelist(n)
fem1 <- aF[el[,1]]
fem2 <- aF[el[,2]]
male1 <- aM[el[,1]]
male2 <- aM[el[,2]]
allFs <- NULL
allMs <- NULL
for (i in 1:length(fem1)){
if (fem1[i] != 0) {allFs[i] <- fem1[i]}
else {allFs[i] <- fem2[i]}
if (male1[i] != 0) {allMs[i] <- male1[i]}
else {allMs[i] <- male2[i]}
}
agesel <- as.data.frame(cbind(allFs, allMs))
sqrtages <- sqrt(agesel)
sqrtages$diff <- sqrtages[,1] - sqrtages[,2]
mean(sqrtages$diff)## [1] -0.23440936 Other Network
6.1 Fit
fit.other <- ergm.ego(egodat_other ~ edges +
nodecov(~ageF) +
nodecov(~ageF<19) +
nodecov(~ageF^2) +
nodecov(~ageM<19) +
nodecov(~ageM) +
nodecov(~ageM^2) +
absdiff(~sqrtage + dat$shift.other*(male==0)) +
nodefactor("deg.mar") +
nodefactor("deg.cohab") +
concurrent(by=~male) +
offset(nodematch("male", diff = FALSE)) +
offset(nodefactor("olderpartnerO")),
offset.coef = c(-Inf, -Inf),
control = x)
#offset(nodefactor("debuted", levels=1)),
#offset.coef = c(-Inf, -Inf),
saveRDS(fit.other, here("fits", "fit.other.rds"))fit.other <- readRDS(here("fits", "fit.other.rds"))
other.netest <- ego.netest(fit.other, DissolutionOther)
summary(fit.other)## Call:
## ergm(formula = ergm.formula, offset.coef = ergm.offset.coef,
## target.stats = m, eval.loglik = FALSE, control = control$ergm.control)
##
## Iterations: 3 out of 200
##
## Monte Carlo MLE Results:
## Estimate Std. Error MCMC % z value
## offset(netsize.adj) -9.903488 0.000000 0 -Inf
## edges -6.065217 1.606071 0 -3.776
## nodecov.ageF 0.448504 0.083211 0 5.390
## nodecov.ageF<19 -0.337883 0.228428 0 -1.479
## nodecov.ageF^2 -0.008739 0.001446 0 -6.045
## nodecov.ageM<19 -0.050555 0.200336 0 -0.252
## nodecov.ageM 0.270125 0.074329 0 3.634
## nodecov.ageM^2 -0.003829 0.001247 0 -3.071
## absdiff.sqrtage+dat$shift.other*(male==0) -3.103887 0.128697 0 -24.118
## nodefactor.deg.mar.1 -4.387958 0.335242 0 -13.089
## nodefactor.deg.cohab.1 -3.673664 0.330056 0 -11.130
## concurrent.male0 -3.374662 0.272150 0 -12.400
## concurrent.male1 -1.471438 0.299298 0 -4.916
## offset(nodematch.male) -Inf 0.000000 0 -Inf
## offset(nodefactor.olderpartnerO.1) -Inf 0.000000 0 -Inf
## Pr(>|z|)
## offset(netsize.adj) < 1e-04 ***
## edges 0.000159 ***
## nodecov.ageF < 1e-04 ***
## nodecov.ageF<19 0.139096
## nodecov.ageF^2 < 1e-04 ***
## nodecov.ageM<19 0.800770
## nodecov.ageM 0.000279 ***
## nodecov.ageM^2 0.002132 **
## absdiff.sqrtage+dat$shift.other*(male==0) < 1e-04 ***
## nodefactor.deg.mar.1 < 1e-04 ***
## nodefactor.deg.cohab.1 < 1e-04 ***
## concurrent.male0 < 1e-04 ***
## concurrent.male1 < 1e-04 ***
## offset(nodematch.male) < 1e-04 ***
## offset(nodefactor.olderpartnerO.1) < 1e-04 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## The following terms are fixed by offset and are not estimated:
## offset(netsize.adj) offset(nodematch.male) offset(nodefactor.olderpartnerO.1)6.2 MCMC Diagnostics / GOF
mcmc.diagnostics(fit.other)
gof.o <- gof(fit.other, GOF = "model")
saveRDS(gof.o, here("fits", "gof.o.rds"))gof.o <- readRDS(here("fits", "gof.o.rds"))
plot(gof.o)
6.3 Static Diagnostic
static.other <- netdx(other.netest, nsims = 500, dynamic = FALSE)##
## Network Diagnostics
## -----------------------
## - Simulating 500 networks
## - Calculating formation statisticssaveRDS(static.other, here("fits", "static.other.rds"))static.other <- readRDS(here("fits", "static.other.rds"))
static.other
plot(static.other)6.4 Dynamic Diagnostic
dynamic.other <- netdx(other.netest, nsims = 5, nsteps = 1000, ncores = np,
set.control.ergm = control.simulate.ergm(MCMC.burnin=1e5,
MCMC.interval=1e5),
set.control.stergm = control.simulate.network(MCMC.burnin.min=1e5,
MCMC.burnin.max=1e5))
saveRDS(dynamic.other, here("fits", "dynamic.other.rds"))dynamic.other <- readRDS(here("fits", "dynamic.other.rds"))
dynamic.other## EpiModel Network Diagnostics
## =======================
## Diagnostic Method: Dynamic
## Simulations: 5
## Time Steps per Sim: 1000
##
## Formation Diagnostics
## -----------------------
## Target Sim Mean Pct Diff
## edges 1951.961 1949.035 -0.150
## nodecov.ageF 48913.199 48823.057 -0.184
## nodecov.ageF<19 2287.271 2285.052 -0.097
## nodecov.ageF^2 1321309.489 1318712.896 -0.197
## nodecov.ageM<19 2181.387 2180.544 -0.039
## nodecov.ageM 52313.794 52230.803 -0.159
## nodecov.ageM^2 1511564.771 1509422.882 -0.142
## absdiff.sqrtage+dat$shift.other*(male==0) 553.355 553.436 0.015
## nodefactor.deg.mar.1 60.275 60.282 0.012
## nodefactor.deg.cohab.1 63.409 66.256 4.489
## concurrent.male0 41.854 40.940 -2.184
## concurrent.male1 146.741 147.834 0.745
## offset(nodematch.male) NA 0.000 NA
## offset(nodefactor.olderpartnerO.1) NA 0.000 NA
## Sim SD
## edges 33.197
## nodecov.ageF 879.840
## nodecov.ageF<19 41.388
## nodecov.ageF^2 26937.054
## nodecov.ageM<19 40.336
## nodecov.ageM 908.185
## nodecov.ageM^2 29079.782
## absdiff.sqrtage+dat$shift.other*(male==0) 13.786
## nodefactor.deg.mar.1 7.517
## nodefactor.deg.cohab.1 7.645
## concurrent.male0 6.279
## concurrent.male1 11.095
## offset(nodematch.male) 0.000
## offset(nodefactor.olderpartnerO.1) 0.000
##
## Dissolution Diagnostics
## -----------------------
## Target Sim Mean Pct Diff Sim SD
## Edge Duration 55.000 51.781 -5.853 51.175
## Pct Edges Diss 0.018 0.018 0.167 0.003plot(dynamic.other)
plot(dynamic.other, type = "duration")
plot(dynamic.other, type = "dissolution")
6.5 Checking Mean Deg Dist by Age/Sex and Age/Sex Mixing are close to data
n <- dynamic.other$nw
# mean deg dist
plot(summary(n ~ nodefactor(~ageM)), type="l", main="Males")
lines(mo, pch=19, col="blue")
legend("topright", legend = c("Empirical", "Model"), col=c("blue", "black"), lwd=1)
plot(summary(n ~ nodefactor(~ageF)), type="l", main="Females")
lines(fo, col="blue")
legend("topright", legend = c("Empirical", "Model"), col=c("blue", "black"), lwd=1)
# age/sex mixing
aF <- n %v% "ageF"
aM <- n %v% "ageM"
el <- as.edgelist(n)
fem1 <- aF[el[,1]]
fem2 <- aF[el[,2]]
male1 <- aM[el[,1]]
male2 <- aM[el[,2]]
allFs <- NULL
allMs <- NULL
for (i in 1:length(fem1)){
if (fem1[i] != 0) {allFs[i] <- fem1[i]}
else {allFs[i] <- fem2[i]}
if (male1[i] != 0) {allMs[i] <- male1[i]}
else {allMs[i] <- male2[i]}
}
agesel <- as.data.frame(cbind(allFs, allMs))
sqrtages <- sqrt(agesel)
sqrtages$diff <- sqrtages[,1] - sqrtages[,2]
mean(sqrtages$diff)## [1] -0.1753593#netests <- list(marcoh.netest, other.netest)
#saveRDS(netests, here("Setup", "FinalNetworks", "objects", "full", "nodefactor", "netests_nodefactor.rds"))